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A single phonon in a nonlinear nanomechanical resonator (NAMR) can block the excitation of a second 
phonon [Phys. Rev. A 82, 032101 (2010)]. This intrinsically quantum effect is called phonon blockade, and 
is an analog of Coulomb blockade and photon blockade. Here we predict tunable multiphonon blockade in 
coupled nonlinear NAMRs, where nonlinearity is induced by two-level systems (TLSs) assuming dispersive 
(far off-resonance) interactions. Specifically, we derive an effective Kerr-type interaction in a hybrid system 
consisting of two nonlinearly-interacting NAMRs coupled to two TLSs and driven by classical fields. The 
interaction between a given NAMR and a TLS is described by a Jaynes-Cummings-like model. We show 
that by properly tuning the frequency of the driving fields one can induce various types of phonon blockade, 
corresponding to the entangled phonon states of either two qubits, qutrit and quartit, or two qudits. Thus, a 
/c-phonon Fock state (with k — 1,2,3) can impede the excitation of more phonons in a given NAMR, which we 
interpret as a /c-phonon blockade (or, equivalently, phonon tunneling). Our results can be explained in terms of 
resonant transitions in the Fock space and via phase-space interference using the s-parametrized Cahill-Glauber 
quasiprobability distributions including the Wigner function. We study the nonclassicality, entanglement, and 
dimensionality of the blockaded phonon states during both dynamics and in the stationary limits. 

PACS numbers: 42.50.Pq, 85.85.+j, 42.50.Dv 


I. INTRODUCTION 

Nanomechanical and optomechanical devices are a versa¬ 
tile technology mo, with a range of applications in the quan¬ 
tum regime M- Examples include small mass or weak-force 
detection IMD, quantum measurements ca, and quantum- 
information processing. The first success in putting a me¬ 
chanical device in a quantum state was performed m by 
purely cryogenic means, due to the frequency of the mechan¬ 
ical phonons M being larger than the thermal energy. Since 
then lower-frequency devices (which thus have larger mass) 
have been put in quantum states using side-band cooling via 
microwave and optical cavities EHEl. These breakthroughs 
have been followed by the observation of state transfer 03 - 
l22l between an electromagnetic cavity and the mechanical 
system, with the goal of developing hybrid mechanical cir¬ 
cuit devices |[T9l[23l[24l for applications in quantum technolo- 
gies MM- 

Moving into the nonlinear regime with such quantum 
nanomechanical devices is desirable for several reasons: It 
will allow us to observe highly nonclassical effects, as well 
as allowing to individually address different transitions within 
the mechanical system, so that it behaves as a mechanical arti¬ 
ficial atom. In the realm of quantum optics IZTtiSbll and circuit 
quantum electrodynamics llT7ti39l , nonlinearities are typically 
associated with “photon blockade” (also referred to as the op¬ 
tical state truncation by quantum nonlinear scissors ll40ll4T]| ). 
In this regime the nonlinear nature of the spectrum of an op¬ 
tical cavity 14^ . induced via, e.g., a Kerr nonlinearity, means 
that the presence of a single photon within the cavity prevents 
the transmission of a second photon. In nonlinear mechanical 


systems one expects an analogous “phonon blockade” ll43ll44ll 
to arise. This requires either strong intrinsic nonlineari¬ 
ties da, or induced nonlinearity via ancillary nonlinear sys¬ 
tems (like qubits or artificial two-level systems 1461 ). Finally, 
while coupling between mechanical phonons and electromag¬ 
netic photons has been achieved in the quantum regime da, 
the controlled coupling between multiple mechanical modes 
has so far been restricted to classical devices dl- Such con¬ 
trollable coupling would enable the observation of entangled 
states dSHi, and the Bell inequality violation with massive 
objects da, as well as the realization of coupled mechanical 
qubits 15^1^ . 

Our goal in this work is to study the combination of 
nonlinearity-induced phonon blockade effects via the cou¬ 
pling between the mechanical modes of NAMRs. We will 
show that the infinite-dimensional mechanical states of the 
NAMRs, under proper resonance conditions, can effectively 
be truncated to the states of finite-dimensional systems of 
either two coupled qubits, a three-level system (called a 
qutrit) coupled to a four-level system (referred to as quartit or 
ququart), or, in general, two coupled d-level systems (qudits). 

While our model can be considered as a quantum limit of 
the classical systems studied in Ref. 1541 . we discuss explicitly 
how the nonlinearities can be tuned using an ancilla two-level 
system (TLS). We will show how this combination of phonon- 
blockade and two-NAMR (or two-mode) coupling leads to 
multiphonon blockade (or phonon tunneling), in analogy to 
the predictions of multiphoton blockade 155H581 and closely 
related photon tunneling l36l[58ti6Ql . For example. Ref. l36l 
provides a pedagogical explanation of how photon blockade 
can lead to the observation of a single-photon tunneling effect 
in an analogous way to how Coulomb blockade can lead to the 
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FIG. 1: (Color online) Schematic diagram for the couplings in the 
dissipative hybrid system described by the Hamiltonian given in 
Eq.Q. The system consists of two linear nanomechanical resonators 
(NAMRs), with frequencies and a pair of two-level systems 

(TLSs), with frequencies driven by classical fields, with fre¬ 
quencies and respectively. Moreover, Qn, /n, Qn, and 
Ji 2 denote the coupling strengths of the depicted subsystems; jn are 
the NAMR decay rates; an (aj^) is the phonon annihilation (creation) 
operator for the nth NAMR, and is the Pauli operator for the nth 
TLS. 


observation of single-electron tunneling. 

We discuss here how the form of the multiphonon block¬ 
ade can be tuned via driving, and verify the resultant highly 
nonclassical states with a variety of measures. 

This paper is organized as follows: In Sec. II, we describe 
a model for a hybrid system of coupled linear NAMRs and 
TLSs. In this section and in Appendix A, we also show how a 
Kerr-type nonlinearity can be induced via NAMR-TLS inter¬ 
actions, and derive an effective Hamiltonian for the coupled 
nonlinear NAMRs. The possibility of observing multiphonon 
blockades in this system is described in Sec. III. We summa¬ 
rize several methods to assess nonclassicality, which we then 
apply in our analysis of phonon blockade in Secs. IV and V. 
We conclude in Sec. VI. 


II. MODEL 


conversion and parametric amplification 1541 . Thus, the total 
microscopic Hamiltonian, representing the system shown in 
Fig.[^ reads (hereafter h = 1 and n = 1, 2): 

H = iTres + + ^JC + ^int5 (1) 

Hres = L] + fn [«n + H.C.] , (2) 

n 

M Q 

Hci = Yj + H-C.] , (3) 

n 

Hjc = y]5n{a„cr++ aj,(7“), (4) 

n 

Hint = Ji2 [«! + H.C.] 

X [a2 exp(iw^^|,i) + H.c.], (5) 

where and a\^ are, respectively, the phonon annihilation 
and creation operators for the nth NAMR, 

an = 

- ipn), (6) 

which are given in terms of the position operator momen- 

in) 

turn operator and frequency ccres of the NAMR. Moreover, 
c^n — \^n){^n\ ~ \9n){dn\ IS the PauH Z Operator for the nth 
TLS, while cr“ = | 5 f„)(e„| (a+ = |e„)( 5 „|) is the quhitlow- 
ering (raising) operator given in terms of the ground (|pn)) and 
excited (le^)) states of the nth TLS; is the TLS frequency 

and (^drv) NAMR (TLS) driving-field frequency. 
The coupling strengths of the subsystems, as shown in Fig.[2 
are denoted by fn^ 9n^ and J 12 . The symbol H.c. denotes 
the Hermitian conjugated term. 

This system, described by Eq. Q, can be realized in vari¬ 
ous ways as a combination of two types of implementations, 
e.g.: (i) the proposal of Ref. 14^ for observing single-mode 
phonon blockade in a driven single NAMR coupled to a su¬ 
perconducting quantum two-level system and (ii) the system 
of two nonlinearly-coupled NAMRs, which was experimen¬ 
tally realized in the NTT experiments (see, e.g., Ref El and 
references therein). It is worth noting that we assumed that 
the interaction iTint is additionally driven at frequencies 

and as given in Eq. (m, which slightly generalizes the 
model applied in Ref. 1541 . We also note that the interaction 
described by Eq. Q conserves the number of excitations, in 
contrast to that described by Eq. 


A. Time-dependent Hamiltonian 

We consider a hybrid system, as schematically depicted in 
Eig. consisting of two interacting driven linear nanome¬ 
chanical resonators (NAMRs), described by the Hamiltonian 
iTres 5 coupled to two driven two-level systems (TLSs, qubits), 
given by the Hamiltonian The interaction Hjc between 
the nth NAMR and nth TLS (for n = 1,2) is described 
by a Jaynes-Cummings-like model under the rotating wave- 
approximation. The interaction iTint between the two NAMRs 
can be interpreted as a combined driven process of frequency 


B. Time-independent Hamiltonian in a rotated dressed-qubit 

basis 


In the following, we assume the following large detunings: 

> yq > 0, (7) 


s a;? s =<> - »!i„ > 0. (8) 

The large detuning Arq^ ^ Pn implies that, e.g., the qubit 
states cannot be flipped by the interaction with the NAMR (see 
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Appendix A). While the large detuning enables us 

to omit, in particular, the terms which do not conserve the 
number of excitations in the TLS-NAMR interaction Hamil¬ 
tonian, as will be explained below. Note that the assumption 
that ccres^ — cssential in our deriva¬ 

tion, and is applied only for simplicity. 

First, we transform the Hamiltonian H, given in Eq. 0. 
into a rotating reference frame by the unitary transformation 


Unit) = jq exp , (9) 

n 

which results in the following effective Hamiltonian 

H' = U^HUn-iUl^^UR, (10) 

which can be given explicitly as 

H' = + + + (11) 

-^res ~ ^ ^ /n(<^n+<^n)5 (12) 

n 

-^q / V 2 2 (12) 

n 

^JC = ^^n[anCr+exp(-iA^^^t)+H.C.], (14) 

n 

^/nt ~ ^ 12(<^1 + <^ i )(<^2 + <^ 2)5 ( 12 ) 


where cr^ = cr+ + cr“. Note that this Hamiltonian is still time 
dependent. 

Now we diagonalize the qubit Hamiltonian, given by 
Eq. ( pA) ), by transforming it into a dressed-qubit basis follow¬ 
ing the method described in, e.g.. Refs. MM- Thus, one 
finds 

= E ^Pn, (16) 

n 

where A^ = y^A^ + and the dressed-qubit operator 
= \En){En\ — \Gn){Gn\ cau be defined by the dressed- 
qubit basis states |[62ll : 

\En) = COS Xn\en) Esin Xn\gn)^ 

\Gn) = -sinXn|en) +COSXn|^n), (17) 

where x^ = (1/2) tdin~^{Qn/An). It is seen that the dressed 
nth qubit refers to the nth TLS dressed with the nth NAMR 
phononic field (for a general discussion see Ref. fSTl ). 

To transform the Hamiltonian Hjq, given by Eq. ( p^ , into 
the dressed-qubit basis, first we note that 
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FIG. 2: (Color online) Schematic diagram for the couplings in the 
dissipative hybrid system described by the effective Hamiltonian 
given in Eq. ( |28| ), which consists of two nonlinear NAMRs driven by 
classical fields. Here Fn and J denote the corresponding coupling 
strengths; Kn is the effective Kerr nonlinearity of the nth nonlinear 
NAMR, and Sn is its energy. Other symbols are the same as in Fig.[^ 


Eq. ([T§, and their Hermitian conjugates, which do not con¬ 
serve the number of excitations. Then, the Hamiltonian Hjq 
can approximately be transformed into 

H'j'c « + H.C.], (19) 

n 

where = gn cos'^{xn)- 

To obtain a time-independent total Hamiltonian, we trans¬ 
form it into a qubit rotating frame by applying the standard 
unitary transformation 

Uq = exp{-iH'^t). (20) 

This results in 

H'j'l: = E9'n {an/o+exp[-i(A(") - A„)i] + H.c.} . 

n 

( 21 ) 

By assuming Arq^ = A^ (for n = 1,2), one obtains the time- 
independent Jaynes-Cummings Hamiltonian in the dressed- 
qubit basis. 


H'^i = Y.g'n{anPl + aip-). ( 22 ) 

n 

Thus, after these transformations, the total Hamiltonian reads 
H"' =HZ + H'Ji + HZ, (23) 

where = H!,, and F"; = Note that F"' = 0. 


C. Effective Hamiltonian with a qubit-induced nonlinearity 

We assume another large detuning 


acr+ = cos^{xn)ap'^-siv?{xn)ap^E\sin{2xn)ap^^, (18) 


^ 9n'i (^4) 


where p^ = |Gn)(^n|, and /)+ = \En){GnV By re¬ 
calling the assumption, given in Eq. ([^, we can write 
acr+ cos^{xn)(ipn • For example, if Vtn/An = 0.1 then 
cos^(xn) = 0.9975, sin^(x^) = 0.0025, and sin(2xn)/2 « 
0.05. Thus, we can omit the second and third terms in 


which leads to dispersive interactions. Note that condi¬ 
tions 0 and ( [24| ) are not contradictory as can be shown as 
follows: By denoting = Vtn/An, Eq. (24) can be given as 


— ^71 


lEr\ 




\Anrl > g'^ 


(25) 









4 


Thus, Eqs. ^ and ( [24| ) can be combined as the following hi¬ 
erarchy of conditions: 


such that 


(26) 


For example, if = 0.1 then we require 0.005A^ ^ g'^. 

Now we describe how the TLS-NAMR interaction can ef¬ 
fectively induce a Kerr-type nonlinearity. To show this, we 
can expand the Hamiltonian H' in a power series of the pa¬ 
rameter 


^ ^ ^ ^ flnCOS^[|tan ^(»n/A„)] 

" VA2+n2 _ A„ ^ 


(27) 


such that I An I ^ 1. As derived in Appendix A, one can keep 
terms of such expansions up to only and assume that both 
TLSs remain in their excited states \En) during the whole sys¬ 
tem evolution, which can be observed for the large detuning 
Arq^ ^ as given in Eq. ([^. Then, the total Hamiltonian 
H'", given in Eq. (23), can be transformed into the following 
effective Hamiltonian: 


ffeff = E [4"el(0, 1) - 4:v4«n 
n 

+J (ai + a|)(a 2 + a^), 
where the Kerr-type Hamiltonian 


(^n ^n) 


(28) 
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FIG. 3: (Color online) Energy levels for the Kerr-type Hamiltonians 
given in (a) Eq. (with n = 1) and (b) Eq. (with n = 2 and 
£■2 = 


where l)£r„, and( 7 ^* = is a constant 

term, which can he ignored. Thus, Eq. ( [^ corresponds to 
1). In the following we will also analyze another 
special case of Eq. ( [32l ) corresponding to 

-^KCTr(l’ 2) = Kn{alan - - 2) -h E^^alan, (33) 


= KnO^alanan + (29) 

describes an effectively nonlinear nth NAMR. Here the nth 
NAMR energy and the effective Kerr nonlinearity are 
given by 


£n — ^res + 2 A’n + 9n^n{^ ~ (30) 

Kn = -9n><l, (31) 

respectively. Moreover, = /^(l + 2^^) is an effective 
driving-field strength, and J = Ji 2 (l + ^A^) is an 

effective coupling between the NAMRs. All these coupling 
coefficients are shown in Fig. A few lowest energy levels 
for the Kerr-type part of this Hamiltonian, given in Eq. (29), 
are shown in Fig.j^a), where we set n = 1. 

Equation ( [28| ) describes an effective driven Kerr-type self¬ 
interaction (an anharmonic model) in one phonon mode non- 
linearly coupled to another phonon mode. In the context of 
various types of photon and phonon blockades, this phonon- 
phonon model can formally be considered as a two-mode gen¬ 
eralization of the single-mode phonon model of Ref. na. 
Moreover, our model with the nonlinear coupling between 
the NAMRs can be interpreted as a generalization of (i) the 
linearly-coupled hybrid model studied in Ref. 1441 to de¬ 
scribe single-phonon and single-photon blockades and (ii) 
the linearly-coupled optical model with single-mode |[63H65l 
and two-mode drivings 163 leading to two-mode single¬ 
photon blockades. 

The Kerr Hamiltonian, given in Eq. ( [2^ can formally be 
rewritten as (/c, / = 0 , 1 ,...) 


l) = K„(aia„ - k)(aia„ - 1) + 

( 32 ) 


where, for simplicity, the term is omitted. A few lowest 
energy levels for this Kerr Hamiltonian are shown in Fig.[^b), 
where n = 2 and £2 = ^ 2 ^. 


The driving of qubits, as given in the Hamiltonian ([^, 
can tune the effective Kerr nonlinearity. In addition, as 
will be discussed in the following and was also observed in 
Refs. 16411651 . there is another mechanism for tuning the Kerr 
nonlinearity in the coupled anharmonic oscillators. Thus, 
even if the NAMR decay rates 7 ^ are much larger than the 
NAMR driving strengths and the latter are much larger 
than the Kerr nonlinearities Kn, strong single-time photon (or 
phonon) antibunching can still be observed as an indicator of 
the photon (phonon) blockade l64l[^ . 


We will analyze free and dissipative evolutions of the 
NAMR systems described by the Hamiltonian given in 
Eq. (^), under two different resonance conditions, as spec¬ 
ified later in Eqs. and ([^). 


III. PHONON BLOCKADE 


A. Phonon blockade in two models 

In this section we will analyze phonon blockade in two 
models, which are special cases of the general model de¬ 
scribed by Feff in Eq. ( [^ under different resonance condi¬ 
tions. 

(n) 

In one model, we assume that the frequency of the 
driving field of the nth NAMR is tuned precisely to the effec¬ 
tive free energy £n (n = 1,2). Then, the effective Hamilto- 
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FIG. 4: (Color online) Energy levels Emim 2 ^ given by 
H\mi,m 2 ) = Emim 2 \mi,rn 2 ), for a simplified model 2 , i.e., the 
system of two uncoupled nonlinear NAMRs, described by the Hamil¬ 
tonian H = ^K^rr’ given in Eq. assuming 8 = 8 i = 82 . 
Here, mi and m 2 are the Eock states ofthe two NAMRs. It is seen 
that the levels E 2 Q and E 02 are off-resonance if Ki, K 2 ^ 812 and, 
thus, they are much less occupied than the other shown levels. This 
explains the occurrence of phonon blockade in this coupled system. 


EIG. 5: (Color online) Energy levels for a single NAMR, as a sim¬ 
plified version of model 2: How to induce a three-phonon resonant 
transition in a single nonlinear NAMR described by the Hamiltonian 
given in Eq. ( [3^ . If the driving field frequency is resonant with the 
transition between the energy levels | 1 ) and \ 2 ), E 2 — Ei — 8 , then 
one can also induce the three-phonon transition between the energy 
levels |0) and |3), since E 3 — Eq = 38. Solid (dashed) lines denote 
real (virtual) energy levels. We note that only single-phonon tran¬ 
sitions can be observed if the driving field frequency is tuned with 
the transition between other levels \n) and |n + 1) (n / 1). More¬ 
over, we cannot observe such multiphonon transitions if the system 
is described by the Hamiltonian given in Eq. ( [2^ for any n. 


nian iTeff. given in Eq. ( [^ , simplifies to 

-^eff ~ ^ — 1) 

n 

+ Fji (<^n + <^n) + J (<^1 + tl|)(a2 + ( 22 ), (34) 

n 

which is referred here to as model 1. The occurrence of single¬ 
phonon blockade in this model is explained in Fig.ffl 

(n) 

In another model, we set the frequency of the driving 
field of the first (second) NAMR to be tuned with the effective 
free energy 8 '{ = 81 {82 = f 2 + 27^2). Thus, the effective 
Hamiltonian iTeff reduces to: 

= Aria|ai(a|ai - 1) + K2{a\a2 - l)(a2<^2 - 2) 

+ Fn {cLn + Cl\^ + J {cLl + Ci\){ci2 + Cii) ^ 


which is referred here to as model 2. Note that we ignored 
in Eq. (35) the irrelevant constant terms (for n = 1, 2). 
Some lowest energy levels for the Kerr-type Hamiltonian of 
Eq. ( [^ forn = 1 are shown in Eig.j^a) in comparison with 
those for n = 2 shown in Eig. [^b). A closer analysis of 
Eig. I^b) shows that E^ — Eq = 3^2 = , which implies 

that a three-phonon resonant transition can be observed, as 
shown in Eig.[^ The occurrence of multiphonon blockade in 
this coupled system can clearly be understood by analyzing 
Fig-IU 


B. Phonon blockade in non-dissipative systems 

In order to show phonon blockade in the non-dissipative 
model 1, we start from the analysis of the system, described 
by the Hamiltonians without dissipation. Hereafter we 
assume that both NAMRs were initially in the ground phonon 
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8 

^'—E^=2K 

FIG. 6 : (Color online) Energy levels for a simplified model 2, i.e., 
for the NAMR system described by the Hamiltonian iTgff» given in 
Eq. ( [^ , assuming two drives and the same Kerr nonlinearities Ki = 
K 2 = K and free energies = ^2 = ^ of both NAMRs, and 
setting J = Fi = F 2 = 0. This graph explains the occurrence 
of multiphonon blockade in this coupled system. The occurrence of 
the resonant three-phonon transition is explained in Fig. We note 
that, in particular, the following levels are off-resonance: Eqi = 
8 7 ^ Eio, Ei 2 — 38 ^ E 21 , Eo 2 — Ell — 28 ^ Eiq + 8 , and 
E '20 = 28 -\- 4iT 7 ^ Eiq + 8 . 


states, = 0)) = |00). A generalization of such solu¬ 
tion for other initial states is simple. For simplicity, we as¬ 
sume that both NAMRs are driven equally with the strength 
El = F 2 = F. Then, under these assumptions, the solu¬ 
tion of the Schrodinger equation for the wave function = 
exp(—is given by 


It/;) = coolOO) + coilOl) + ciollO) + cii|ll), (36) 
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with the time-dependent probability amplitudes: 

Coo = ^ g4iFt ^ , 

Coi = CIO =--exp(-iJt)sin(2Ft), 

cii = Coo - exp(iJt), (37) 

as calculated for simplicity in the two-qubit Hilbert space. 
Only for short evolution times, these solutions approximate 
well our precise numerical solutions, which were obtained in 
a high-dimensional Hilbert space and plotted in the left frames 
of Fig.|7] Much better agreement with these precise numerical 
solutions can be found by diagonalizing the Hamiltonian 77'^ 
in a two-qutrit Hilbert space. Unfortunately, we cannot obtain 
a compact-form analytical solution in this case, as discussed 
in Appendix B. 

As a measure of the quality of phonon blockade (or phonon 
truncation), one can calculate the fidelity, defined as F{t) = 
Emi,m 2 =o.i 1 ^. where the probabilities are 

computed precisely in a large-dimensional Hilbert space. For 
the same parameters as in the left frames of Fig. we find 
that the fidelity periodically oscillates between the values 
0.977 and 1. This shows that the evolution of phonons in the 
NAMRs is practically confined in a two-qubit Hilbert space in 
model 1 even without dissipation. 

Let us now analyze phonon blockade in the non-dissipative 
model 2. For simplicity, we again assume that the driving 
field strengths are the same, Fi = F 2 = F. The solu¬ 
tion of the Schrodinger equation for the wave function I'lp) = 
reads, 

= coo|00) + co 3|03) + cio|10) 

+ ci3|13)+C2i|21)+C22|22), (38) 

where the probability amplitudes can be found only numeri¬ 
cally, as explained in Appendix B. We plotted the evolution of 
these probability amplitudes for model 2 in Fig. (solid blue 
curves) in analogy to those shown in Fig.[7]for model 1. 

Analogously to model 1, we can quantify the qual¬ 
ity of phonon blockade in model 2 by calculating the fi- 
delity F{t) = E(mi,m 2 ) kmi,m 2 p; but now {mi,m 2 ) = 
(0, 0), (0, 3), (1, 0), (1, 3), (2,1), (2, 2). For the dissipation- 
free evolution shown in Fig.[^ we find that the fidelity F{t) G 
[0.9643,1]. Thus, we can conclude that the evolution of 
phonons in the NAMRs according to model 2 even without 
damping is effectively confined in the Hilbert space of an en¬ 
tangled qutrit-quartit system. 


C. Phonon blockade in dissipative systems 

In the standard description of dissipation under Markov’s 
approximation, the evolution of the reduced density operator 



Jt Jt 


FIG. 7: (Color online) Non-stationary phonon blockade in model 1: 
Evolutions of probabilities P{mi,m 2 ) = |cmi,m 2 l^ for the non- 
dissipative (left frames, 7 = 0 ) and dissipative (right frames, 7 = 
J/3) systems of the two nonlinear NAMRs described by the effective 
Hamiltonian 77'^^ plotted for the rescaled time Jt, where J is the 
interaction strength between the NAMRs. Probabilities P(?tii, m 2 ) 
for the other values of mi, m 2 are negligible on the scale of these 
figures and, thus, not presented here. We set the Kerr nonlinearities 
Kn — lOJ, the drive strengths Fn — J, and the mean number 

— in') 

of thermal phonons = 0.01 (right frames), for n = 1,2. It 
is seen in the right frames that the oscillations are rapidly damped. 
Surprisingly, the stationary damped states are highly nonclassical as 
it is shown in other figures. 

p{t) is governed by the master equation, 

P = -i[He«,P\ 

+ (Saipan - anttip - /)a„4) 

n 

+ ^ ^ 3“ — a\^ajip — paj^a^), (39) 

n 

where 7 ^ is the nth NAMR decay rate (damping constant), 
^th ^ = {exp [cc/ ( /cbT)] — 1 } “ ^ is the mean number of thermal 
phonons interacting with the nth NAMR, ks is the Boltzmann 
constant, and T is the reservoir temperature at thermal equilib¬ 
rium. For simplicity, we assume equal decay rates 71 = 72 = 
7 and mean thermal-phonon numbers = nth- In 

our numerical analysis, we focus on the steady-state solutions 
Pss = p(^ ^ 00 ) of the master equation, obtained for p = 0 . 
We obtain such numerical solutions by applying the inverse 
power method implemented in Ref. 1^ . 

Examples of dissipative evolutions of phonon-number 
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FIG. 8 : (Color online) Non-stationary phonon blockade in model 2: 
Probabilities P{mi , m 2 ) same as in Fig.[7]but for the non-dissipative 
(blue solid curves, 7 = 0 ) and dissipative (red dashed curves, 7 = 
J/3) NAMR systems described by the effective Hamiltonian 


probabilities are shown by the red curves in Fig.[7]for model 1 
and the red dashed curves in Fig. for model 2. It is seen 
that short-time oscillations are rapidly damped. The resulting 
steady-state phonon-number probabilities are shown in Figs.|^ 
and[^ respectively. In panels (a) of these figures, we plot¬ 
ted the probabilities P(mi, m 2 ) = (mi, m 2 |pss|^i: ^ 2 ) for 
the two-NAMR density matrices pss- While in panels (b) and 
(c), we presented the single-NAMR phonon-number proba¬ 
bilities = {m\p^^\m) for = Tr 3 _^(pss), with 

n = 1,2. 

We can clearly interpret these results as single-, two-, and 
three-phonon blockades corresponding to the cases shown in 
Figs.[^b),p^b), and[^c), respectively. 


We note that the steady states are not pure, contrary to the 
standard assumptions made in analogous studies of single¬ 
photon blockades of optical 165][67l and optomechanical ISSll 
systems. Indeed for the examples of the states shown in 
Figs. 0 and [T^ we found that their purities are the follow- 
ing: = 0.4212 and TxipH^f = Tripif)^ = 0.5670 

for model 1, and Trpgg = 0.1471, Tr(/9ss^)^ = 0.3920, and 


Tv{pis^)^ = 0.3212 for model 2. 

Finally we note that, for simplicity, we applied here the 
standard master equation, given by Eq. assuming two 
separable dissipation channels for the NAMRs. A more pre¬ 
cise description, which could be especially important for a 
stronger coupling between the NAMRs, should be based on 
a generalized master equation within the general formalism of 



(b) 




photon number m 


FIG. 9: (Color online) Stationary solutions for model 1 de¬ 
scribing an effective two-qubit system and, thus, corresponding 
to single-phonon blockade: The phonon-number probabilities (a) 
P{mi,m 2 ) = (mi, 77121 Pss I mi, m 2 ) for two NAMRs and (b) 
P^^\m) = (m|Tr 2 (pss)|m) for the first NAMR (and, equiva¬ 
lently, (m) for the second NAMR) for the steady-state solutions 
Pss of the master equation ( [^ with the Hamiltonian H'^^, given 
by Eq. ( [^ , assuming the same parameters as in Fig. [7] More¬ 
over, (xy) denotes all phonon numbers such that x,y > 1 , so 
P{xy) = 1 - Paucl (b) clearly shows 

the occurrence of single-phonon blockade in every NAMR. 


Breuer and Petruzzione (see sect. 3.3 in Ref. (691). In this 
approach both NAMRs dissipate into usually-entangled dissi¬ 
pation channels. An explicit form of such generalized master 
equation for two strongly-coupled infinitely-dimensional sys¬ 
tems will be presented elsewhere Col. Note that a generalized 
master equation for an infinitely-dimensional system strongly 
coupled to a qubit system has already been well studied (HI. 


IV. PHONON BLOCKADE IN PHASE SPACE 

Here we apply the Wigner function W = and the 

Cahill-Glauber s-parametrized quasiprobability distribution 
(QPD) in order to visualize the nonclassical properties 
of the phonon-blockaded states studied in Sec. III. 

The Wigner function for a two-mode (or two-NAMR) state 
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(a) 

0.2 



photon numbers (WpOTj) 


(b) 

¥ 


Q, 



photon number m 


(C) 



photon number m 


FIG. 10: (Color online) Stationary solutions for model 2 de¬ 
scribing an effective qutrit-quartit system and, thus, correspond¬ 
ing to multiphonon blockade: The phonon-number probabilities (a) 
P(rni^m 2 ) = (mi,m 2 |pss|'^^i,'^^ 2 ) for two NAMRs, and (b,c) 
— (m|Tr 3 _n(pss)|m) for the nth NAMR (n = 1, 2) for 
the steady-state solutions pss of the master equation ( [^ with the 
Hamiltonian i^eff ^ given by Eq. ( [^ , and the same parameters as in 
Fig. and (xy) denotes all phonon numbers such that x > 2 and 
y > 3. Panel (b) [(c)] demonstrates the occurrence of two-phonon 
(three-phonon) blockade in the first (second) NAMR. 


p can be given by 

FFi 2 (q^ 1 ,0^2) = ^12(^1,Pi, ^2,^2) 

= {<11-^1^Q2-X2\p\qi+Xi,q2+X2) 

X exp [ 2 i{piXi +^22:2)] dxidx2, ( 40 ) 

in terms of the canonical position and momentum op¬ 
erators, and an = Qn P Wn for the nth NAMR. It is seen 


that Eq. ( [40| ) is a straightforward generalization of the Wigner 
function for a single-mode (in our case single-NAMR) case, 

w{a„) =W{qn,Pn) = ~j Xn\pn\qn + Xn) 

X exp {2ipnXn) dXn, (41) 

where pn = p^'^^ can correspond, e.g., to Ti^-nP for n = 1, 2. 
Specifically, the single-NAMR Wigner function W{an) can 
be considered as the marginal functions of the two-NAMR 
Wigner function W (o^i, 0 ^ 2 )- 

Figure pH] shows the Wigner function W ( 0 : 1 ) = W ( 0 : 2 ) for 

the steady state pss = Tis-npss iP = 1,2) for some cho¬ 
sen values of the coupling and damping parameters. Unfortu¬ 
nately, the contribution of the vacuum state in pss is dominant, 
as shown in Fig.|^ Thus, the Wigner function for looks 
like a slightly deformed Gaussian representing the vacuum. 

To show this deformation more clearly we also plotted the 
s-parametrized Cahill-Glauber QPD, W^^\an) for 5 = 1/2. 
For simplicity, we analyze this QPD only for a single-mode 
(i.e., single-NAMR) case, while the extension for the two- 
and multi-mode cases is straightforward. The Cahill-Glauber 
QPD W^^\an) can be defined in the Fock-state representa¬ 
tion of an arbitrary-dimensional single-mode state p as fol¬ 
lows 1721 : 


= y] {k\p\l){lp\an)\k), (42) 

/c,Z =0 

where 

(43) 

for s e [-1,1], c = :lexp[-2|Q!„|2/(l - s)], Xa„ = 
4|a„|V_(l - s2), y = 2/(1 - s), ^ = (s + l)/(s - 1 ), 
and I/f * are the associate Laguerre polynomials. As for the 
Wigner function, an is a complex number, where its real and 
imaginary parts can be interpreted as canonical position and 
momentum, respectively. The operator T^^\an) is defined 
in the Fock representation by Eq. ( [43] ). In the special cases 
of 5 = —1,0,1, the QPD W^^\an) becomes the Husimi Q, 
Wigner W, and Glauber-Sudarshan P functions, respectively. 

Here, we apply a well known definition of nonclassicality 
(see, e.g., Ref. ca and references therein): A photonic or 
phononic state can be considered nonclassical if and only if its 
Glauber-Sudarshan P function is not positive (semi)definite, 
which means that it is not a classical probability density. Thus 
only coherent states and their statistical mixtures are classical. 

The Wigner functions and 1/2-parametrized QDPs for the 
steady states of models 1 and 2 are shown in Figs. im an dEH 
respectively. These steady states pss are nonclassical, as they 
correspond to the partially-incoherent finite superpositions of 
phonon Fock states, which are not mixtures of coherent states 
(in particular, they are not the vacuum). The nonclassicality 
of these states pss is clearly seen in the non-positive functions 
pjy(s=i/2) negative regions are plotted in blue). How¬ 
ever, the nonclassicality of pss is difficult to deduce from the 
plots of the nonnegative Wigner functions. 
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Re(a^) Re(a^) 


FIG. 11: (Color online) Stationary solutions for model 1 as in Fig.|^ 
but for (a) the single-NAMR Wigner function W{ai ) = W ( 0 ^ 2 ) and 
(b) the single-NAMR quasiprobability distribution (QPD) function 
— W^^\a 2 ) with parameter s = 1/2 for the steady-state 

(n) 

solutions pss = Trs-nPss (for n = 1,2). Note that the negative 
regions of the QPD functions are marked in blue. 


To understand this apparent discrepancy, we recall a known 
relation between two single-mode QPDs, and for 

the chosen parameters s < sq CH: 

=C' J exp 

(44) 

where c' = 2/[7r(so — s)]. This relation means that any QPD 
for s < So can be obtained from by mixing it with 

the Gaussian noise. In particular, the Wigner function can be 
obtained in this way from the P = and QPDs. 

By decreasing the parameter s from sq = 1, the QPD 
for a given nonclassical state, becomes less and less negative, 
and finally becomes nonnegative at some s' > —1. As a re¬ 
sult, in the analyzed examples shown in Figs. [TT] and the 
negativity of the P-function for pss is only partially lost in 
the QPDs but completely lost in the corresponding 

Wigner functions. 


V. ENTANGLEMENT, DIMENSIONALITY, AND 
NONCLASSICALITY OF NAMRS 

To analyze more deeply the nonclassical properties of the 
generated phonon states in the two NAMRs in models 1 and 2, 
we apply the following measures of quantum correlations: the 
negativity and its closely related entanglement dimensionality, 
as well as the entanglement potential, as a measure of nonclas- 
sicality. 


A. Entanglement 

To quantify the entanglement of a bipartite state p of arbi¬ 
trary finite dimensions, we apply the negativity N, which can 
be expressed as 1741 


N{p) 


ll/lli-l 


(45) 



Re(a2) Re(a2) 


FIG. 12: (Color online) Stationary solutions for model 2 as in Fig.[T0l 
but for the single-NAMR Wigner functions W{an) (a,c) and the 
single-NAMR QPD functions W^^\an) with parameter s=l/2 (b,d) 
for the steady-state solutions pis^ = Trs-nPss (for n = 1,2). 


via the trace norm 111 1 of the partially-transposed statistical 
operator p^. This entanglement measure is closely related to 
the Peres-Horodecki criterion. The negativity N is an entan¬ 
glement monotone and, thus, can be used in quantifying en¬ 
tanglement in bipartite systems. However, the negativity does 
not detect bound entanglement (i.e., nondistillable entangle¬ 
ment) in systems more complicated than two qubits or qubit- 
qutrit M. The negativity can be interpreted operationally. 
For example, the logarithmic negativity. 


-^cost (p) = log2[(V(p) + l], (46) 


quantifies the entanglement cost under operations preserving 
the positivity of the partial transpose (PPT), which is, for 
short, referred to as the PPT entanglement cost 17511761 . 

The evolutions of this entanglement measure are plotted in 
Fig.[^a) for model 1 and Figs.[^a) for model 2, by includ¬ 
ing and excluding the dissipation. The oscillations of F^cost (p) 
are rapidly damped; however the entanglement is not com¬ 
pletely lost in the infinite-time limit. Indeed, for the coupling 
parameters AT, J, F, decay rate 7 , and thermal-phonon mean 
numbers nth specified in the figures, the entanglement be¬ 
tween the NAMRs is found to be Acost(Pss) = 0.1413 for 
model 1 and almost three times smaller Acost(Pss) = 0.0494 
for model 2. 


2 
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B. Dimensionality 

The negativity also determines the dimensionality I^ent of 
entanglement, which is the number of degrees of freedom of 
two entangled subsystems. Specifically, the entanglement di¬ 
mensionality T)ent ior a bipartite state p is simply related to 
the negativity N{p) as follows fTl\ 

Dent(p) = 27V(/9) + l = ||/||i. (47) 

More precisely, the least integer > gives a lower bound 
to the number of entangled dimensions between the entangled 
subsystems of p fTTl . According to Eq. ( [47] ), = 1 for 

separable states {N = 0). This measure could be useful for 
characterizing even a single test system (in our case, a single 
phonon mode) with unknown quantum dimension. This can 
be done in a standard way “by entangling [the test system] 
with an auxiliary system of known dimension and measuring 
the negativity, a lower bound to the number of quantum lev¬ 
els in the test system can be found” fTTl . In our case of two 
NAMRs, we can directly apply the negativity, without the use 
of an auxiliary system, to determine a lower bound to the num¬ 
ber of quantum levels in the total system (see also Ref. (781). 

The evolutions of the entanglement dimensionality are plot¬ 
ted in Fig. p^b) for model 1 and Fig. p^b) for model 2. 
Since the entanglement dimensionality and the entanglement 
cost are closely related, we can conclude, the same as for 
Ecost{p)^ that T)ent does vanish in the steady states. Specifi¬ 
cally, the entanglement dimensionality between the NAMRs 
reads: T)ent(Pss) = 1.2058 for model 1 and T)ent(Pss) = 
1.0696 for model 2. 


C. Nonclassicality 

The negativity can also be used in quantifying the nonclas¬ 
sicality of a single-mode photonic or phononic state pn via 
the so-called entanglement potential (EP), which is defined 
as |7l|79l 

EP{pn) = log 2 {A^[exp(—jiJi)(/9„(8)|0){0|)exp(iiJi)]+l}. 

(48) 

Here H = + a„6’^) describes a balanced beam splitter 

or a linear coupler, where and b are the annihilation opera¬ 
tors of the input modes. The basic idea behind this measure in 
optics is as follows: If a single-mode nonclassical (classical) 
photonic state is combined with the vacuum at a beam split¬ 
ter then the output state is entangled (separable), for which 
various entangled measures (including the negativity) can be 
applied. By generalizing this concept for phonons it is enough 
to interpret this ancilla beam splitter as a linear coupler. 

The evolutions of the nonclassicality of single NAMRs are 
plotted in Fig.[^c) for model 1 and Figs.[^c) and[^d) for 
model 2. We find the following nonzero values of EP(pss^) 
in the corresponding steady states: EP(pss^) = EP(ps?) = 
0.1126 for model 1, while EP(pss^) = 0.1354 for the first 
NAMR and EP(pss^) = 0.1770 for the second NAMR in 
model 2. 



Jt 


FIG. 13: (Color online) Non-stationary solutions for model 1: (a) 
Entanglement, measured by the PPT entanglement cost Ecost(p), 
(b) dimensionality of entanglement 71dim(p), and (c) nonclassical¬ 
ity, measured by the entanglement potential EP(p^^^), of the first 
(and, equivalently, second) NAMR for the states p generated in the 
non-dissipative (blue solid upper curves, 7 = 0 ) and dissipative (red 
dashed lower curves, 7 = J/3) systems described by the effective 
Hamiltonian i^eff • Parameters are same as in Fig.|^ 

Here, for brevity, we studied only the evolution of one non¬ 
classicality measure. In future work, it might be physically 
interesting to compare it with the evolution of other mea¬ 
sures ca and witnesses (SOI of nonclassicality. 


VI. DISCUSSION AND CONCLUSIONS 

We studied tunable phonon blockade, which can be intu¬ 
itively understood as follows: Any number of phonons can 
be generated in a harmonic resonator. However, this is not 
possible in an anharmonic resonator, which is characterized 
by nonlinear (nonequidistant) energy levels. So, if the driv- 
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Jt 


FIG. 14: (Color online) Non-stationary solutions for model 2: Same 
as in Fig. 


13 


but for the effective Hamiltonian i/eff- Additionally, 
panel (d) shows the nonclassicality, measured by the entanglement 
potential EP(y9^^^), of the second NAMR. 


ing field is in resonance with the transition between the two 
lowest levels (say |0) and |1)), then it is not in resonance with 
the transitions between the other levels. Thus, single-phonon 
blockade can be observed. We showed in detail that a higher- 
order n-phonon blockade can also be observed in a dissipative 
nonlinear system if the driving field is resonant with the tran¬ 
sition between other levels |0) and |n). By applying coupled 


nonlinear systems, instead of a single system, one can more 
easily tune various types of multi-phonon two-mode block¬ 
ades, as studied in detail in this paper. 

It is important to clarify the main differences between 
our model of coupled oscillators and that studied in, e.g., 
Refs, im |64j [65l: (1) Here we assumed that the oscilla¬ 
tors are nonlinearly coupled as described by the Hamiltonian 


^int ^ (til + tii)(ti 2 + <^ 2 ) contrast to the linear cou¬ 
pling given by i^int ^ (a|a 2 -b aiaj), which was applied in 
Refs. EUlMlEa. Note that this linear coupler (which is for¬ 
mally equivalent to a beam splitter or a frequency converter), 
does not change the nonclassicality of a total system Ca. In 
contrast to this, the nonlinear coupler can increase the non¬ 
classicality of a total system, as measured by, e.g., the entan¬ 
glement potential. It is also worth noting that, by applying 
our precise numerical calculations, we found the steady-states 
of the NAMRs to be only partially coherent (partially mixed), 
while the steady states calculated in, e.g., Ref. ll65]| were as¬ 
sumed to be completely coherent (perfectly pure). (2) We 
have derived an effective Kerr-type Hamiltonian from a mi¬ 
croscopic one, given in Eq. ([B, while an analogous Kerr-type 


Hamiltonian in Refs. Il44ll64tl65]| was assumed without deriva¬ 
tion. Moreover, we studied here the Kerr interaction under 
two different resonance conditions, as described by Eqs. ( [29] ) 
and(^). References naiMiisi discussed only the Kerr in¬ 
teraction given by Eq. ( |29| ). (3) We analyzed the blockade 
of mechanical phonons, contrary to photon blockade studied 
in Refs, m isa. Moreover, the interplay between single¬ 
phonon blockade in one oscillator and single-photon block¬ 
ade in another oscillator was studied in Ref. HH. Here we 
predicted tunable /c-phonon blockades (with k = 1,2,3) in 
each oscillator, where the /c-phonon Eock state impedes the 
excitation of more phonons. To our knowledge, multiphonon 
blockade has not been studied before. 

In conclusion, we showed here a rich tapestry of phonon 
blockade effects in two coupled nonlinear nanomechanical 
resonators. Different types of phonon blockade could be 
“picked out” of this tapestry by controlling the nonlinearity 
via ancilla TLS, and by changing the driving frequency of 
the resonators themselves. Within these different types of 
phonon blockade, the coupled NAMRs can be made to be¬ 
have like two coupled qubits, a qutrit coupled to a quartit, or 
even two coupled qudits. We verified this picture by looking 
at the nonclassical properties of these states including their 
single-NAMR nonclassicality, and two-NAMR entanglement 
and entanglement dimensionality. The nonclassical properties 
of these states were also analyzed in phase space by applying 
the s-parametrized Cahill-Glauber quasiprobability distribu¬ 
tions and, in particular, the Wigner function. We expect that, 
if realized in experiment, the ability to operate in these dif¬ 
ferent regimes will have a range of applications in quantum 
information and quantum technologies. 
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Appendix A: Derivation of the effective Hamiltonian 


Here we show how to derive the effective Hamiltonian 
given by Eq. (28), from the Hamiltonian H'", given in 
Eq. ([23]). The latter can be divided into the following parts: 


-Ll „ „ “T -O “T 


(Al) 


where 


u'" — u'" 

^sys — ^0 


H 


JC — ^ + Hjq, (A2) 


^drv — 


(A3) 


Here is given by Eq. (22) and 777'^ = is given by 
Eq. 

Our derivation is based on the method described in Ref. IMI 
for the exact diagonalization of the Jaynes-Cummings model 
with the following unitary transformation 


Un = exp[-A(A„)(aJ,/5„ - a„p+)] 


(A4) 


where the operator A (An) = — arctan(2An\/^)/(2\/]Vn) 
is given in terms of the total number of excitations Nn = 
ahcin + \En){En\ in the nth NAMR and TLS. Thus, the 
expansions of the annihilation operators, = U}^aUn and 
Pn — ^nPn^n of R^f. lIMIl Can be rewritten for our dressed 
qubit states as follows: 


(^n — ^in(l + \KiPn) + Pn + ^n^nPn + C^(^n):(A5) 

Pn = Vn + KanPn “ + C>(A^), (A 6 ) 


where = 1 — k)^(a\an + 1/2) and O(A^) denotes the 
omitted terms of order ^ A^ and higher. Now one can easily 
transform the Hamiltonian H'" into 


H = uIuIh"'U-2,Ux. (A7) 

In particular, by applying Eq. ( |A5| i, transforms into 

fn[an + a^]- (A8) 

n 


If the qubits remain in the excited dressed-qubit states \En), 
given in Eq. ( [Tt] ), then 

{E^E2\Hd„\E^E2)=Y,Uan + ai){l + \l/2) + 0{Xi). 

n 

(A9) 

The assumption of a “frozen” state of both qubits is physically 
justified for the large detuning Arq^ » fln, as specified in 
Eq. (|^. Moreover, transforms into 

^int = ^ 12(^1 + a|)(a2 + ^ 2 ), (AlO) 


and, thus. 


{E^E^\H,^,\E^E^) = J{a^ + a\){a^+a\) + 0{\t,). (All) 

where J = Ji 2 (l + Ai/2)(1 + A 2 / 2 ). Analogously, by gener¬ 
alizing the results of Ref. ID for the dressed-qubit operators 
one can find that 


ffsys = (A12) 

= ^ 0 ' - 5 - v/l+4A^nA2)p- 

n 

~ 77q + ^ ^ ‘^KnCLnCLn + hnPn 3“ -^n(^n) ^nPn^ 


where Hq' is defined in Eq. (A2), hn = Xn(<^n<^n + with 
Xn = ^n^n(l — A^), and the effective Kerr nonlinearity reads 
En = -g'n^n = -{d'nY/ Moreover, U = U 1 U 2 and 
= 0. Thus, one can write 


(7^l7^2|77sys|^1^2) — '^^Cn DndliCin 

n 

+Kr,{aifal + 0(A4), (A13) 

where Cn = Xn/2 (which can be neglected as a constant 

term) and Dn = £n - = A„ + 2A'„ + Xn- Finally, 

the effective Hamiltonian 


77eff = (^1^2|(77sys + ^drv + Tfjnt) 17^17^2 ) , (A14) 

is given explicitly by Eq. ( [2^ , where the terms r\j 0{Xi) are 
omitted. 


Appendix B: Probability amplitudes in Eqs. ( |36) and ( |38) 

The probability amplitudes Cxy{t) = {xy\'ip{t)) (for x^y = 
0,..., 3), given in Eqs. ( [3^ and ( [38] ), can be obtained using the 
eigenvalue decompositions 77'^17^^) = 7^^|7^^), as 

c-xy{t) = '^eyi.v{-iE'J){E'J)0){xy\E'J, (Bl) 

n 

and analogously for 77''f^|7^^') = 7^^'|7^^'). To simplify these 
problems, let us limit the dimension of the Hilbert space to 
that of two qutrits and assume 7fi = 7^2 = 10 J and J = 
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= F 2 (as in Figs. [7] an d[8j. Then the Hamiltonians, given 
in Eqs. ( [34| ) and ( [35) ), reduce to 
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(B2) 


(B3) 


respectively. It is seen that the matrices, given by Eqs. ( |B2[ ) 
and ( |B3[ ), differ only in their diagonal terms. Unfortunately, 
even in these special cases, it is very unlikely that exact analyt¬ 
ical compact-form solutions of these eigenvalue problems can 
be found, as they require finding the roots of sixth and ninth 
order equations, respectively. Especially, highly-irregular os¬ 
cillations of Cxy{t) for model 2, as shown in Fig.[^ confirm 
this conclusion. 
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